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1 ghMEROruCTIGN 


skuméen (1972) claims that progress in numerical modeling 
of the general circulation has been to some degree dictated 
Se wie dest by the rate of development in the field of 
comouter tecmno | oa. Beowever, the limited ability to 
parareterize the effects of small-scale processes in terms 
or large scale motions has been an é@qually important 
wemcting, Pector. Fssentially, the maior problem of numerical 
me@eling of the general circulation is simply that of 
procucing a very long range numerical weather forecast. 

Certainly the equations used in the models must de more 
sophisticated to include those physical processes-which -are 
Ipeemepertent for e short range forecast, bsbut may become 
Cec i@l was the lemeth of the foré€cast is extended. Another 
area where concentrated efforts nave improved the forecast 
invelves the computational techniques employed v0 
approximate and solve the governing equations of the models. 

The Motivation tehind this thesis is to investigate the 
meio Of da m@elatively new computational techniaue to 
Mee Weclc of numerical weather prediction. The finite 
elerent method, long established in engineering, has been 
seriously considered only during the cast deéecede in 
meteorology. This method has great potential for application 


Speetmesogeric preciction models. 





A. BACKGROUND 

Phemmost commen aurerical integration procedure for 
B@aener orediction has been the finite difference method in 
werem the derivatives in the differential equations of 
mettem ere repleced by finite difference approximations at a 
@eiserete set or peiats in space and tire. The resulting set 
Gameomatiors., with appropriate restrictions, can then be 
solwea@ by aleebraic methods. Until récently, the finite 
difference method has been the workhorse in atmospheric 
preeciction models, from their first computer implementation 
mo Lhe present. 

With the introduction of each new Pewerat1 on Oe 
cemputers. the Zap between mUnerical <ferecasts aad 
dtmosomeric observations has decreased. The rate at which 
this gen decreased has slowed down and appears to be 
Wemelbing Off. This would indicate that computer technology 
acon 22 Sune pra@mary obstruction to setter aurerical 
ferecasts. In fact, blegger and faster computers alone kave 
@emonstfated their inability to significantly improve the 
pumerical forecast. 

Rag werterTple, a major limiting mactor Ont, fea 
@ameecerce aporoximations is the truncation error. The 
Nem lomeak Weather Service 7 Layer Primitive :quation Model 
(7LPE Model), operational from 1966 to 1980, had truncation 
errcers which increased at a rate proportional to tne square 


Om the "“erid spacing. That is, the smaller the grid interval, 


ie 








eee sgeeiier the’ trtmcation error. To increase its accuracy 
mere Weare increasing the grid ratrix density. This would 
mewerre increase@e computer steerage and comrputaticnal time. 
Beewe Or the art computers are capable of vreovidine these 
edé@itional resources. 

The problem now goes beyond numerical techniques and 
Somputer technology. Operationally, the National weather 
service is not cavable (due to monetary restrictiozs) of 
vrovicing a censer concentration oz atmospheric 
Oemertidtous. Therefore, with the osresent density cf initial 
data (observations) and objective analysis techniques 
(eettise the data for grid poirts by interpclating from 
observed data sources), reducing the grid spacing further oa 
the 7LP& Model does not Significantly intrease the accuracy 
cH them solution. 

This additional computer capability can not de utilized 
using finite difference methods. Therefore, new numerical 
integretion techniques must be investigated, such that siven 
tae sem densitre of observed data, superior sciutions are 
produced. 

Two alternative techniques, the spectral method and the 
fimase el@rent rethod, have started to gain attention. foth 
the spectral and finite element methods require more 
computetionel time per forecast time step than dees the 
finite difference method. For example, the finite element 


Menmoa Meaquires dn @€quation sclver to invert a lerger matrix 


1S 








ee “ech tie step for each variable. In this sense, these 
methods were held back by computer technology, but recent 
Somwomces i= computer teckhtolozgy (i.e. larger and faster 
storage cevices, multi-processors, ets} have made these 
aivemnative numerical techniques competitive. 

tor longs range weather predictions, the spectral method 
éovolied over the globe or hemisvnere is a natural methoc, 
Ge te the existence of efficient transforrs for the 
aonlinear terms on spherical geometry. It also eliminates 
Uwe trwm@eetion error for the horizontal space derivatives 
and the nonlinear instability (aliasiag). For these reasons, 
zlobal spectral models have been developed and implemented 
On an operational level, replacing the global finite 
Gitference models. 

Fowever, because the spéctral harmonics are glotally 
Bowmer thaqg locelly defined, it is thought that for problems 
Co “meme e@eteiled limited area forecasting, the finite 
elerent method is more suitable. Pioneering work to adapt 
finite element methods to meteorological epplicaticns has 
Giiea’ dome by Cullea (1973,1¢974 and 18979), Staniforth and 
Mitchell (1977), Binsman (1975) and Xelley (1976). The rost 
peg@mt f¥aite element meteorological rodel at the Naval 
Postgeredwate School was written by Kelley (197€) with the 
Ceumeeommension of Er. R.v. Williams. It is this study that 
Vor esermumcs @ basis for this thesis. The model written ty 


Beivey wili be altered and used for comparative testing with 








improvec finite element forms implemented by this author. 
some of the techniques and codes developed by Kelley are 
ee ec SmploOved in this thesis. Older (1981) develored a 
Becoiniaque to smoothly vary the grid geometry in the domain. 
This technique is also implemented both on <Xelley’s model 
oe Wet) )«6©6the|6Cnew6OCU formuletion to give greater versatility 


weer testing the model perfcrmance. 


BS. C3JECTIVES 
Tmeeoojectives for this thesis can be divided into twe 
GeummeTries: i) meteorolosy, 2) computer science. First, the 
meteorological objectives of developing imprcved finite 
elerent forms for shallow water equations are as follows: 
Peo (Cllttere (1922) afer collaboration with Dr. 
M.J.F. Cullen, showed how equilaterally shaped elements 
pre@rcecd sienificantly Gemner “sesuilts “than did ‘other 
triangular elements. ZYelley (197€) used right triangular 
elem@nts in the implemertation of a two-dimensional finite 
element model using the primitive form of the shallow water 
equations. 4 considerable amount of small-scale noise was 
GOuserveca In the solution. EFereafter, this model, which was 
develoved vy Zelley (1976), will be referred to as the 
promi ti ve rodel. This frrs¢ objective involves 
re-implémenting the orimitive model using equilaterally 
shaped elements and comparing the results to those in 


Yelley’s thesis. 





2) - ‘Silliams and Zienxiewicz (1981) presented new 
rare @lement tecnoniques for formulations for the shallow 
Memes Cua taens, whict use differently shaped “unctions to 
eppreximrate the different ceéevendent variables, which ina 
effect Stagger the variables. Schoenstadt (1980) 


deronstrated the advantage of spatial Staggering of 


fd. 


ependent variables Poet te «6cit ference models. The 
emoircation of this technique to finite element rodels is a 
mereret ‘Sxtension. and excellent results were obtained by 
Williems and Zienkiewicz (1951) from application of these 
mem formulations on linearized one dimensional cases. The 
omyecetave here ts to implement the new forms on the 
mreeitive rodel and azaina cdo quanitative comparisons of the 
results. 
Zz) - The maitor emphasis in this study deals with the 
Bo t@rer etionm aad comparison of the vorticity divergence 
pore 0: tite’ shallow water equetions, which is descrited in 
oeet) i Shapter 21. This formulation hes the following 
Gawenteees. First, the geostrophic adjustment process is 
Peetee Yetter than in the primitive form of the eauation;ns. 
meeondly, the velocity and height fields are evaluated at 
ee seme erid poiit, where the pesteprimitive form recuires 
Stageering these dependent variables. And thirdly, a larger 
me st@p Ws allowed dwe to the semi implicit form of the 
equation. Again comparisons between the results from the 


Weumecity ciwerzence aad vrimitiwe model are presented. 





-me computer science aspect of this thesis was primarily 
cevoted to the implementation of the different models and 
Bee Segte ane architecture of the program. Finite element 
momeemes TedQuire more computationel time than do finite 
efferexce femeges, DOwe only im™=the solution of the 
GEeuettons. Sut also in the amount of corputation required to 
eee uate each term in the equations. 

"Se Miplememtations of tnese two dimensional models, 
elthough complex when viewed from the surface, have a lot oS 
generelity ame “Mecumeamicy in the operations required. 
Yersetile modules can be written to ease the implementation 
eae fecilitate changes. The objective here is to efficiently 
implement these new forms and demonstrate the utility of 


teese versatile modules for future implementations. 


c. -2a 1S STRUCTURE 

Ris thesis presents the results obtained fror tests of 
the variovs finite element formulations. The results are 
compared to those from the primitive modél written by “Yelley 
iecme). Becompanying the results is a detailed discussion of 


cme eetormitetion and implementaticn process. 
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CHommerm [1 of tme thesis presents a tutoriel of tk 


(D 


finite elerent metnod and the area coordinates system used 
Pe one Seelwetion of the elerent integration. The Galerkin 
finite element method used in all the rodels is developed 


eee cemaiec bo the advection equation in one dimension. 
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G) 


ot 
fay) 


eeweer Iii presents the detailed description cf the 
verticity-divergence shallcw water model. Zere the ecuations 
are = cn Sec Wee Reece te «€©6Galerkin tethnocd. A 
mo seussiem of the computational technique used ic presented 


along with the rodel’s paysical varameters. 


emepter IVY cresents a desc 


rs 


iptive overview of the 


SomMPuter implementation. The 


QO 


hapter includes a list of 
Sorvems emeilable for testing, a brief description of the 
mame? Comoacticn techaniaue and the formulations of the 
mersatile modubes used to implement the complex ecuations. 


Cha 


‘oJ 


ters V through VII discuss the results obtained fror 
the different experiments. Chapter V briefly describes the 
Deeomeuive model usred for all comparisons and the results 


“rom cwanging the element shape to equilateral triangles. 


Chenter VI reformulates the primitive model so that the 


OY 


Mov0tGmtial is staggered with respect to the velocity 
Seueenyum@. “Or simplicity, the continuity equation is also 
Seeearizec. Chacter VII comvares the results from the 
vorticity-divergence mcdel developed in Chapter III to those 
PoP the primitive model. 

The lest chapter summarizes the results from all the 
Ewpetamests anc identifies what areas reauire follow oa 
wort. The source code for the vorticity-diverzgence model is 


wee Sented in Appendix A. 
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Meets GO: ten the case with an original develonprent, it is 
momen CGierricult to quote an exact date on which the frinite 
eeeement fethod was invented, but the roots of the rethod can 
Me tracec tack to these separate zroups: suena 
MaeMematicians, paysicists and eéngineers. Since the early 
€eveloprenats of the finite element method, a large amount of 
research has been devoted to the technique. However. the 
finite element method obtained its real impetus by the 
in€ependent cevelopments carried out by engineers. I[ts 
Seseateat SamMplicity in toth ohysical interpretation and 
mathematical form has undoubtedly been aS much- behind its 
Mopuleri1ty as is the digital computer which today permits a 
Meatistic solution of even the most comolex situations. 

The name ” finite element “was coined in a vater dy 
2.W. Clough, in which the technique was presented for plane 
stress analysis, as discussed in Eathe (197€). While finite 
element methods have mace a ceen impact via the field of 
Olle) m@ehnamics, wnere’ it can ve said that today they 
Pesrecemt the gexnerally accented method of discretizing 
comtimuur croblers for computer-based soluticn, the same 
aeeee@c-06€© SOUmtlhCUCtCOlCUCte SClCdtrue in fluid mechanics or atmospheric 


Deediction. 





Numerous finite element formulations are CuUErel Gly 
Seriatle. Strang (1973), wNerrie (1973) and Zienkiewicz 
P4271? present cetailed theoretical discussions of each. The 
Galerkin method. tke most popular finite element rethod, is 
cesecribed a Petal velow “eha> used in “the equation 


* 


ermuletion later. 


he indt 


r=4 
ta) 


LEMENT CONC 


taj 


pS 

“we problem of solving partial differential equations 
came ce specified im one@ of two ways. In the first, finite 
difference methods specify the depnendcent variables at 
eertein grid points in space and time, and the derivatives 
are evaluated using Taylor series approximations. Secondly, 
ume celeulus of wvwariation requires the minimization of a 
functional over a domain, where a functional is defined as ae 
variational integrél over the domain. 

Peeecatcoulus of Variation apvroacn creates a purely 
eaFsicél mwrodel where the functional equivalent to. the known 
OPE Terentieél eGquetions are Kaown. Its major disadvantage is 
ee) i) «(Cla tS the method only to tBose problems for which 
Peeclomals ezist. Finite element methods, an extention of 
fees “Mathod. derive mathematical anproximaticns directly 
fror the differential eauations goveraing the problem. The 
@eWentage here is that it extends the method to a range of 
mror ems for which a fumctional may not exist, or has aot 


beer discovered. 


20 





Phe finite clement method <¢ivides the domain inte 
Se@ememns Of finite elements (usually of the same form). 
mO7esS are iocated along the boundary of the elements, 
Me@eity at the element vertices and at strategic positions 
amrese@e, centroid, etc.) in the interior and on the sides 
of faces of elements. 

Commonly usec elements are triangular, polygonal or 
polyhedral in ferm for two-dimensional problems. The choice 
G@- @iements ceperes of the type of problem, the number of 
elerents desired, the accuracy required and the available 
comvuting time. To begin with, the element must be able to 
pPepreseat derivatives of wp to the order required in the 
Semmeteom procedure, and to guarantee contiauous mirsit 
cerivatives aero sis Ehe element Houmcaries, to avoid 
Singulerities. Trianguler elements are employed in this 
Testis wvecause they can be used effectively to represent 
irregular boundaries, aad/cr geometry, and also to 
Seopeentrate coordindve functions in those regions_of-_the 
Gomein where repidily varying solutions are anticipatec. 

Consider the problem o? solving approximately the 


differential equation 
Llu) = fix) II-1 


“abe re Deieoemee Geer fienemitiel ovierator, u the dependent 
Maneewlo ane £(x) i6 a specified forcing function. Suppose 


feet =f) (UGS tolCUCbe:6SOclved in the domain a = x = b and that 





Smerooriate boundary conditions are provided. The residual R 


is formed “ror [I-11 as follows: 
oie tie) =) > lIl-z 


@remetitiical step is to select a trial farily of 
approximate solutions (the members of a trial family are 
Ormen celled tasis functions). The basis function S 
prescribed ‘functionally) over the domain in a piecewise 
fashion, element by element, and are generally a combination 
of low order ovoolynominals. A one dimensional example is 
shown in Figure 1, wherein the domain (x axis) is divided 
into six elements (line segments) A through F. The basis 
Reectfons dre linear and one 15 shown for node + only in 
Figure 1. The function has a value of unity over node 4, and 
Cecreases linearly to is zero at nodes 3 and 5 and Zero 
elsewhere. 

Come der a series of linearly independent basis 
eemcti@m@s VY: -(%), as in- Figure 1. Now u(x) .can be 


3 
aeporoximavec With @ finite series as follews: 


where ?, is the coefficient of the jth basis function and has 
d Value equal to u at mode j. 
Suestitutime this approximate solution II-3 wherever u 


dppears in the differential equation II-i 


_ & 


L(V, ) acer) a= 3 II-4 








*uoTpouNT 


b 


Sy ceo eee eel ra 


me 


' «i » 


"T amndty 





she best solution will be one which in some sense 
merece se trae residual R to a minimum value at all points in 
moe comain. =y defiritios, the residual obtained using the 
exact dirferential equation is identically zero everywhere. 
The residnal &, formed in equation [I1-4, is minirizied when 
Mubtipliec with a weighting function, inteagrated over the 
doméin and set equal to zero. This process is knceown as the 


weighted residual method 


{ RW ax = 3 tI- 


cn 


where W is the weighting function and is referred to as the 
test function in the following development. The weighted 
residual method minimizes the errors of the residuals, such 
that the summation of all the positive and negative errors 
ede te zero. 

The Galerkin method, the most povular finite-element 
method, is more general in application and is a special case 
of the method of weighted residuals, as discussed by Pinder 
and Gray (1977). The requirement imposed on the weighted 
residual method forming the salerkin rethod is: 

> twes test (weiehtime) function be equal to the 
basic (trdal) fumctiow ¥ = V . This process 
leads in general to the best approxrimatiioa of 


tien so lution. 
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The final Galerkin form is obtained Dye souosti@eting [1-4 


into II-§, yielding 


» D 
{ H2edyr 23 a [ estas = 2 ies 
a a 


pif this procedure is reveated for N points in the dorain 


ad system with N equations and N unknowns will be generated. 


3. GALEPKIN APPLICATION 
me Tollowing example taken in vart from Haltiner and 
Williams (1$€@) applies the Galerkin rethod to the advection 
equetion with linear elements 
du du 
—-+c—=2 11-7 
ot ox 
Pages @0@ation is dependent in both time and space. The 
treatment of time variation aS impostemt” «for rost 
Mereoroalofical prediction problems. The Galerkin rethod is 
mem copriea to the time dependence because it is fmrore 
Soewewrent to use fiaite differences in time,as is done with 
aes exemple later. The same treatment is applied tc the 
proenostic eauations later, where two finite differencing 
methods are employed to do the time integration. 
Tne Salerkin procedure represents the cependent variable 


eet )C«Cw CC SIM «Clo Ff )60CUtUnCtions that have the prescribed 





Seerial structure as in Figure 1. Approrimate u(x,t) with 


Ree faite series as follows 


wiget) = 


i PA =< 
ane 


j=1 


Beere the coefficient B Cn mcsaumetion of tite, is tre 
Seomeae Value of u at anode j. The basis functions, V(x), are 
mumcetioms Of Spece only and j equals i te 7 for the exarple 
eeetaeure i. The repeated subscript in this forr implies a 
cur ever the repeated subscript. 

the Galerkin equation for the advection equation II-7 is 


octained by setting L = c(d()/dx) and sutstituting in the 


aeoroemim@ate solution [1-8 wherever u is found. 


db 
NX ag. N OV. 
2 —! J vytyex + ¢ § p. aes = 2 II-9 
j=1 Ot j=l om 
a a 
moere i = i to ON, Vs Poem nest. functionmward V, the basis 


function. The domain of integration is given by 4 =x © b, 
eee ee Botegration is doze ia a piecewise fashion, elerent 
dy element. 

In this one-dirmensionel case, @n Equaticn like Ii-¢ is 
writter for each node, i. Considering node 4, what are the 
Messi tle non—Zzero contributions from equation II~9? Figure 
Eeeintewravece the basis and the test function interaction 


awrine the piecewise integretion process. from the 


00 
QO) 





Basis Function 


figure 2. Basis and 
during the piecewise 


Test Function 





test function interaction 


integration process. 


Zt 





Sreeeetron ct the basis and the test function, locally 
Betetee aS wiity at no@e j and linearly decreasing to zero 
es i 6 6and)6=6zero elsewhere, the OnLy Noa 2e 50 
contributions are made when ‘ = 3 over element C, { = 4 over 
elements © and T, and } = & over element D. 


mae CMaluationme of II-S form i = m. which is given in 


Ealtiner and dilliams (1981), leads to the equation: 


ick (a + 4u + u eo em b: -u ) = 2 II-12 
at m m-1 aes m+ 1 m-1 

The bSoundary points, which in this example are nodes 1! 
and 7, are evaluated in the same way as the interior nodes, 


oem the @rception that cyclic conditions are imposed. 


ene time discretization of Il-l@-is done using a ~fi'tite 
ere sertence scheme. Applying leapfrog time cirferencing gives 


ae following G@cuation 








n-l n=1 Beeeeena ) antl — nei 
say Unt 1” aS m m=-1 “n-1 ) 
eae, 26-4. ) = 2 id 
1 m 
Zip ei 


Ime Tesultant equation set, in matrix forr, contains an 
New Meterr where N is the number of nodes. 

The tipfaimis 1 t i of2 tT reir. one-dimension to two its 
Mathematically identical. The domain is now subddivided into 


finite areas, which are triangles in this implementation and 
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the besis functions are linear. However, now they are 
pyrarid skared with value unity at the center and decrease 
t9 zero et tne surrounding nodes, and are zero elsewhere. 
mmeemee c SmowsS this basis function for node 28 outlined in 
heavy black. The value at any node again can be approximated 
by II-3, where j ranges over all nodes connected to node i 


mekverme +: itself. “he connectivity for node i = 28 in 


tay 


feure © is j = 15,16€,27,28,2¢S,59 and 42. 

The integration is still over the entire domain. With 
meen (th basis and the test function zero over the domain, 
except locally over each element, the global integration can 
be performed ty integrating locally over each element. By 
@efinition. this integration can be expressed as an inner 


preoouet of both functions (i.e. basis, test) as follows: 


A 
Usimg this definition and the repeated Suse ri pt 


notétion equation II-S becomes 


pening 1c M464. .V > = 9 11-13 


Peere the dot implies differentiation with respect to time, 
and ioe “Second subscript implies differentiation with 
respect to the second subscript. The local integration may 
Pe Calculated @irectly from exact expressions derived from 


erea coordinates described in detail in the next subdsection. 


N) 
©) 





[\{\/\/ 


—- 


p< = 








Figure 3. Basis function for node 28. The 


shaded area is the complete basis function 
and the V5 , where j = 15, 16,27,28,29,39,40 
are jth node basis functions for node 28. The 


dashed line at node 28 has length unity. 


30 





In summary, the Galerkin procedure involves Sube@ividixrze 
the dorein into finite elements, approximating the dependent 
Wervebies by a linear corbination of lew order ovolynomials 
anc substituting them into the equetions. The ecuation is 
rultinvlied ty a test function, integrated over the entire 
Gomein anc finally the resulting system of equations is 


Selvec. 


C. AREA COORDINATES 

While the Cartesian coordinate system is the natural 
G@eeice of coorcinates for most two dimensional problems, it 
is not convenient wnen working with triangularly snaped 
elements. I[t is therefore necessary to define a special set 
of normelized coordinetes for a triangle. Area, or natural 
Coordinates as they are commonly called, reduce the 
fermiageble task of integrating products between the basis 
mac test fumctions amd their derivatives over a triangular 
elerent and result in easily computable and exact 
Szrpressioas. 

The following development is taken in part fror the 
Towmmietion by Zienkiewics (1971!. Consider the triangular 
element illustrated in Figure 4. There is a one-to-one 
Correspondence between the Cartesian coordinates (X,Y) and 
the area coordinates (Li ,U5,h, ) for the element. Let A 
denote the area of the triangle and Ay A,and A, the areas of 


the subtriangles in Figure 4 such that aA = Ay tA, *A3. 


ol 





(Oxon 1) 


y (X55¥5) 
3 
P a... > ax) 
(G70) 
P(X,Y) or 
G,) 
(1,0,0) 


Fig 4. Cartesian vs. area coordinates 





Fig. 5. Transformation to area coordinates 
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res Geratioushipn between a point P(X,Y) in Cartesian 
coordinates and P(L,, La, L,) in area coordinates can be 


Geen Se the following transformations 


rsx + 1% * Lk 
Y = L,Y + iY + L,Y II-14 
Seat a bon me 
where L, = _ , L5 _=2 aad L3 pes 
A A A 
anc 
es eee O74 2 dt) 72a 
ce (2A + box + aot) /24 aa 


+ 


L. = (2A + bak a3Y) Ven 


myere 2A is twice the area of the triangle and the a’s and 
b’s are defined as in Figure 5. 
i Psweworth motimge that ewery tuple (1 


L L 


ee 2’ 3 
Semresponas to a unigue pair (X,Y) of Cartesian coordinates. 
To Fievre ¢, Lae 1 at vertex i and @ at vertices 2 and @. A 
linear relation exists between the area and Cartesian 
coordinates which implies that values for L,vary linearly 
over the triangle with a value one at vertex i and a value 
@® wero at vertices 2 and 3} and similarly for Lane L.. 
This demonstrates how each component in the tuple (L,, L., L3) 
behaves over the triangle as do the linear basis and test 


functions over the element, as was seen in Fizure 4. Clearly 
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where 5 ns) a wlanear funetiion of the Cartesian coordinates 


[aee. basis, test). 
Zienkiewicz (1971) shows that it is SOsci Le uo 


Pategrete any ovoolynormial in area coordinates using the 


sirole relationshin 


mon.pD rin! Pp! 
Ly Lo ay = CCA Pee 
3 (mn +n + eee) 


maere ti, 2 and p are positive integers and A is the 
elementary area. For an example of this integration 
Deemnsaqwe using imager product aotation, equation [I-le is 


evéluated as follows 


; c! 3! @! A 
Vi axdy = ———————-———_ 2A = - i= j 
(2 +3+ G+ 2)! E 
4 
Cr. A> «= 
_- 1! 1! @! A 
ce ee Ok eras 
J Cue)! 12 
A 
Tite 


The differential operations in area coordinates follow 


ieee Gly “rom the cdi*ferentiation of (II-1&) where 


3 5 by 0 
=i = or II-1¢9 
Ox {=1 2A OL, 
and 
3 5 ay 0 
ey i=l 2A OL, 





as explained earlier (see EFquation II-16), ee ee ear 
function ‘i.e. basis, test) which equals a component L, of 


Oae 2ree coordinate tuple. Therefore 


ov, 6 if i # j 
— = [I-21 
OL, eee et. = 
Conseauently OV, Ox fore = lois 
mo a) lf OM, do dV; b3 OF, by ; 
"3x = — = = 51 =a oe ~ cael I{1-22 
x A Ly 2A L, 2A L. ZA 
0 g 


mS <n ©xaMple, consider the inner product ieee? at 


Mmervwces j} = <, 1 = 1. This integration is evaluated as 


7% 
AY 3 WV. axcy 
2x if Oh 1 


A oes, 


N 

<3 
«<? 
\4 

ib 


b I 
2A (1 +8+ 2+ 2)! 


2 


Nie 


See@refore any inger product in the formulation can be 
Beedgily e@veiuated uSime ar@a coorcinates. Another benefit of 
Men@e this coordinate system is that all of the inner 
procucts are functions of space only and need be computed 


@nly omce. 
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fil. SHABEGN WATER MODEL 
*he gcverning equations fer this model are derived by 
Memes Several simplifying assumptions on the primitive 
\egudtions of motion, which then give the barctropic shallow 
meter €QWatiens. fowever, as mentioned previously, tne 
Shaliow water equations describe many significant features 
of the large-scale rotion of the atmosphere, and therefore 
heve teen usec ia numerous experiments over the years. 

The verticity-divergence form of the equations has 
several advantages. Williams (1981) has shown that the 
meest®oOomic adiustméent process is treated much oetter with 
tee vorticity Givergence formulation than with a direct 
meeetment of the primitive form or the shallow water 
eGuations, such as was used by Kelley (197€). This 
formulation also allows the. velocity components and the 
hé€ight to te carried et the sare nodal points, whereas the 


meet Seobeme for the primitive form of the equations requires 


in 


tapgering of the fields, as seen in Schoenstadt (1980). The 


Worticitv civergence “ornm of the equations is also 
convenient for tne apm ca ti070 of Semi hice t 


@ifferencing, which saves considerable computer time. 


oE 








BR. GOVERNING ECUATIONS 
eye primitive form of the shallow water equations in 


Mew] sicr coordcinetes is 


O74 3 Q 
— +Dd = - —(fJu) - —(fy) an 
Ot Ox oy 
du ag Ox 
a oe ~— + OF - — ee se 
ot ox Ox 
Ov og OK 
a 8 eo 
dt dy Ov 


Bauetiom ‘I1I-1) 1s the continuity equation and the III-2 
ape IilI-3 are the momentum equations, respectively. The 


war iabples are cefined as follows: 


x,y -— the spatial coordinates of the domain 


Uyye— cCompon@ats of the wind vector 


t Bopoveitial = (gravity x free surface height) 
? —- mean geopotential = 49,09¢ meters’ /seconds< 
t time 


IN 


- kipeticvrenerzy 


€@ - absolute vorticity = (3 + ?,) 

8 relativeevorticity 

ee eeriolis .faree (mid-channel ?-plaae) 
D ed vergeace 


The snallow water equations can be written in 


Worticity divergence form as follows: 
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— + Db = - —(pu) - —(pv) III-4 
z p (pv 
es 3 
—— = - —(10) - —(ve) TII-5 
ot Ox ey 
aD 3 3 
— +f = - —(ve) - —(uc) - ve III-6 
ot Ox Oy 


meee CU LLU i =4 is 6Uthe seme continuity ecuation as III-1, III-& 
Ms the vorticity equetion and III-€ is the divergence 
equation. 

s3ecause of the Wervicity 2iverzeice form: of the 
equations, it becomes necessary to solve the time dependent 
Bersepres S and 2 in terms of YW, the strear function 
Ometational part of the wind), and X, the velocity potential 
Seamergent cert of the wind). The initial fields for the 
model will be in terms of WY, %& and #. 

The follewizg diagnostic relationships are defined and 


meee ilat@r in the solution of the equation set. 





u = i XY Ce 
a X* Ky III-8 
where the subscript imolies differentiation 
Coe 2 
Lae 
k= kinetic energy, ks 
Zz 
“ = west <)>, Liao 
ve = v(S + &, IITI-11 
a= wu, il eee 
B= fv, aoe 





Se 77 ¥, Pela 


= Sh fila 


m ACUATZON FORMULATION 


+3 


n 


iD 


GelerKir Method described in Chapter [II is now 


}+- 
Ou 


aemlie to eauations II1i-4 through III-1&. For ease of 
comprehension, the shorthand inner product notation es in 
meete Will be used to simplify the €quations. The detailed 
Galerkin formulation will be shown for equation III-7, the u 
eomeonemt of motion. The rethod follows directly from the 
meerole in Chapter [I of this thesis, which in turn follows 
in vart from Kelley (1976) and Ealtiner and Williams (1981). 

Consider eqvation III-7 and assume that each variable u, 


Yand XT is avvroximated by 


oe, oe. 

"3 5 

y= wie = 
y ¥, 7, mee 


meere the repeéetéed sutscripts indicate summation over the 
mee 06O”™6CUCUHE)6C SSCript. Substituting these approximate 


feemnetons into III-7 yields 


3 
u.V¥. = —heoy. ) oe — hey .) TII~17 
oD dy J 9 ax J J 
Senoe Only the vasis function : MoncomaItecr yond of —SDace, 


eae, ey Se PUrGher simpWified by factoring out the tire 


mt coefficients. 


fu 
m 

cS 
‘ay 
t3 
fd 
(D 
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Pre ext Step requires multiplying by a test function V5 


mo wetscussed in Chapter II, and integrating over the erea 


ff if = 
a Nea = a : on { 
us aVy aA ¥ Vs cA 


oy 


Comain 


A A 
av. 
een Se da roe 
J Ope 
A 


mere See! «form in inner product notation is 


site: > ie o Ue »V5> te OS a eee Nez aes 


7 


Meere the double subscript implies differentiating with 
meemeect to the secomd subscript. 

me three prognostic ecuations (II1I-4, III-& and ILI-é) 
Emewcimilarly advanced using the Galerkin technique to 


wecome, respectively: 


<g.V. i ee = — “ony. .v.> - <B.Y. eerie 
B.V.,V> + ) <I iV 54 V5 OV VG — SBSV aye Vy lS 24 


Gaya ane> = —<(uG) 7, ,¥.> = <{vO).7..,V, DIII-21 
apes oo pee, ok cya 
“~ , = wer ‘ = 7 _ ¢ : \ t Vv > 
<T.V,,9,> <2, ant <(vQ),V,. 5Vy> AUG) Ta Wy, 
+ £.9°7.,V.> ores 
3 S| ge 


meee w is the Leplacian operator and the dot impiiss 
differentiation with respect to the time decrendence in 
meee, TlI-S and ItI~€. 

semi larly, Galerkin equations are formulated for 


smeerctoms L1li-7 through IIfI-1is. 
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oe SGits DISCRETIZATION 

preeeewetion set [1I-2@, III-21 and III-22 is arranged 
fo oat all tne terms on the left nazd side can be treated 
Meolicitiv, and all the terms on the right hand side can be 
mmeaeed explicitly. The explicit time integration will be 
meme cy the leapfrog dirference method. To start the time 
Wee@eenetion, two forward nalf stéovs are taxen, after which 
mee full leastrog scheme is used for the remainder of the 
morecest veriod. 

"he vorticity equation III-21 1s solved indevendenatly 
meer 6CUlL I —2SlhCUm «= CT ITT -22. However, ieiek— 20) and Piet ee 
(Continuity and divergence equations, respectively) are 
Comeumec. TO Explicitly solve etther, decoupling of the 
Bametions is necessary. In this thesis this is doze througn 
wee Oreaic suvstitution of III-22 (solved for D{n+1})) into 
eee o. Once the tire integration is verformed on III-2@, 
P-22 cam be solved “or D(nt+i) using the 4 (nti) value. 


nme #inea! wrnediction equations are 


n+1 P 
I reg SV 5yVgy> + CKV5.¥5>] = 


gy’ 
1 
enol = (srry) ° 
? ae - 
a ce x »V 5? =a Aus RY 5%, 
- $i ie, elie? + Ce 
Z _ _ Nn 1 
FA aes) 5 <v Va V < (ul) 3<W Vy ms 
- 2; oa ee ey 
OD 
aa a) G ¢ ane age 
3[ ( buys F Vaso 4 > + (DV) s<Veve V5 > 
III-23 





MaeresA = 4/{20t) , B = A/d , C = B/(2at) ard [BDRY}] 1s the 


Peeeorurooiic Coundery contribution, see Section 3. 


n+i n-1 


$5 °<V5-7> = 35° <¥a.¥,> 
~ Bet [lus recht > + (VOR Kw Ty >] 111-24 
= “Y..7,> = 1, eae + (at/2) (Or ket y Ae 
aS fees j et 
BO kv Wo, Ve > + 2(vQ)aCVa Vg > 
- 2(uQ)s <Va,V 4 > + Eee Vs Vy >) III-25 


iter these three elliptic equations are colved, the 


tre 


meereory Of the variables [II-7 through III-15 is updated. 

A lerge time step can be applied to this form of the 
@eeeloW water ecuations due to the semi-implicit nature of 
bee Gauetions. “his is very important since finite element 
methocs senerally require more computer time per time steno. 
meme, Vorticity-divergence forrulatica acts as a filter, which 
Slows Ccwn the high frequency waves in the soluticn. The 
twe-Gitensional advective stability criterion for a linear 
element. derived by Cullen (1973), was used to determine the 
morrect tite step, 


4X 


Icl Vé 


Ni 
& A 
ta 
bt 

l 
nw 
Dv 


4t 


rere st is the time step in seconds, 4x the shortest gric 


svecing in meters and c the fastest phase velocity. 





ee COMPUTATIONAL TECHNICUSS 

Phe “fiteal vrognostic equation set requires the solution 
of @ Felmholtz esuaticn for 0 and Poisson equations for Y 
and 7. The most common method of solution used by 
meteorologists has been the successive over relexration 
method (SOR) in which an tnitial guess of the solution is 
mecde and then progressively improved until an acceptable 


fevel of accuracy is reached. SOR is employed in the 


solution of the equations,where II1I-23 can be represented ov 


71 ee =a ey 
and III 24, III 258 by 
seid) = {0} II I-28 


where 7 tie baplkacian operator, {M] = Se Tatrix, {x} 
= tme dependent variable in vector naotation, © - constant as 
in III-22 and fb} the right hand side of the equation or the 
Percing function. 

The mass matrix (M}], dimensioned (axn), is a matrix of 
coefficients whose rows are the equations of the system to 
ve solvec. There exists a one to one sorrespondence tetween 
Pmewrows 6F the mass watrix and the nodes of the domain. 
Beach equation has a term (colurn) for each node, where a 
Ger 8eto term represents connectivity. todes are connected 
if they are both vertices of the same element. Obviously [M] 


Poa eeperse Fatrix conteining the inner products for the 


£5 





left hand side. Chanter 4 of this thesis will describe the 
Metrix compaction vrocedure. 

The forcing function {bt}, dimensioned (Sele iivel ves 
Only variables at the current time step and is easily 
computed using four very versatile subroutines descrited in 
@ereil iszithe next chapter. 

pre tmitial guess to start SOR is the orevious time sten 
solution. Aa average of 308 passes per equation are needed 
mon @@e> time step. The solution is consitered to have 
converged to its #inal value when the residual for each node 
has been reduced to some acceptably small value. 

Wee Giagrostic equations III-7 through III-15 must also 
be solved every time sten. Fowever, the same techniare is 
me used for these eouations. Dr. “.J.P. Cullen suggested an 
WMreer Yelaxetion scheme for which three passes over the 
@omain should produce a solution of acceptable accuracy, 
Bemom vie coefficient matriz is so strongly ciagonally 
Serinent. “ass lumping ef the coefficient. matrix is used for 
the first guess. This technique requires revlacing the mass 
metrix ““] by the identity matrix ‘I]. A first zuess of this 
Teme is 2ble to describe most of the large scale features, 
meron am tttirn rieduces the numoer of iterative passes over 
the field. Successive passes converge to solutions which 
describe smaller scale motion, apvroximately to the same 
order of magnitude as introduced dy computational errcr, $0 


tagt further iterations are not neeced. 





B. GRID GEOMETRY 

Bee Scemaiz of this model is a cylindrical channel, with 
meemes Lemeth of £945 Km and wicth of Z8@S Km. The channel 
memeee tes a belt around the earth and itt vroves to be an 
Geoellent test bed for comparing with the finite element 
formulétions used by Xelley (1S7€) and Older (1981). 

ewe domain is suvdivided into equilateral triersles as 
meow in @ignre €. Most of the test runs for this thesis use 
@eeeric mesh Which has 15€ nodes and 288 elements. This 
Deeeecmentetion is not restricted to one grid vattern. The 
techniaue developed by Older (1681) to vary the nodal 
georetry smoothly to achieve areas of denser and coarser 
Peoeimtton issalso implemented, as in a third grid pattern 
eee woerie’s the nodal geometry abruptiy. A short discussion 
c- these nodal geometries with accompaning iliustratiors of 
mayen is oréesented in Chapter VII, where the different test 
cases are described. 

Grovicecontiauwhty 25 «dssmred in the x @irecticn oy 
wravvinge the domain around the earth to forma cylindrical 
@emaan. this has the advantage of eliminating the east-west 
Peemcertes anc it simulates the flow aroundc the earth. The 
Semey COuncdries cn this dofrain are the norta-south walls and 


~ 


emeir treatment will be discussed shortly. 
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tet Tre, CONDITIONS 

Ser (™eationec previously, the reformulation of the 
meverning Seme tloaus into tae Wor be loyacalrerceence 
Shallow water equation set requires solving the tire 
cepencent variadles in terms of the stream function and 
Pemeccity potential. The continuity equation is not altered, 
so that its solution is expressed in terms of #. 

For the basic testing of the models performance, simrvle 
waewemetaic «€©6SinuSOidal initial conditions are used to insure 
mmeerost accurate analysis possible and to simplify the 
Comparisons. 

"he Sinusoidal initial fields are graphically shown in 
Figure 7 &@s Z-dirensional surfaces. The geopotential field 2 
Gemsist of a half sine wave’ in the y directicn and a sirgle 
mesame weve th the x direction. The stream function Y, 
Seaeculapec by dividing the geovotential field by fae 
coriolis force, has the same physical structure as %. The 
meunccers DoteataAal % Bas Sa single .sinre Wave -in the x 
Gerw@ction and a Half sine waive in the y direction. 


hese 1n7%ial coM#ditions are computed as follows 


= to ASina, cosa 


Z ml 
3 = Of, Pitas 
xX = Csina,sine, Cuaiseeeeostroonic civergence 
where A - arbitrary amplitude 
i? We coriolis velue for mMid—-channe!l latituce 
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3-cGimersional view of the inital fields. 


a) @ and Y initial fields. 
ingeial field. 
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mean flow 
Ya te ce latitude value of y 
= emeene tree =ceonpoteitial heieht 
= 46,200 m@/s * 
= ny /a 
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ro - Wave 1urter 

W Channel width 

EL - channel length 

Co = =(2j,0axyBa)/ (25 + $5) 
3 0 af - os 


G. sCUNDARY CONDITIONS 

Boumeary conditions are only required on the north and 
meoutme wells of the grid domain. Due to cyclic continuity, 
the domain is wrapped around creating a cylinder eliminating 
Bee est end west boundaries. However, careful attention to 
@e@%ell is n@edéed during the implementation to assure this 
continuity. Separate boundary conditions are applied to each 
ee. She pRecictor equations Il1-23, I[1i-24 and III-25. These 
Gom@itions aré comtuted for the wall nodes only and are 
eeplieée @iriag each pass through the relaxation scheme. 

The vorticity eauation III-24, the most sensitive of the 
peee@ietor Soldtions to colve, requires Y on the north-south 
boundaries to remain constant for the entire forecast 
ceriod. ‘Since this equation is solved in terrs of Y, the 


fritial rorth-south Y values are saved and assigned to the 


49 





Mmeume@erY petuts efter each pass through the relaxation 
Sweroutine. 
mie Oreper boundary condition for the divergence 
Semetion I11-25 would be OX/dn = 2. Eowever, for the purpose 
muneeeewies Stay, there is more interest in the sinusoidal 
merretion in the y direction and not In the region of the 
malic. smerefore X = 2 is approvriate. 
me Comtinui ty Semcon si ll=23, the most complex 
mmeurectoOr equation, requires that there be 10 mass flux 
through the north-south walls. The geéeostrophic ooundary 
Reme ition 
az 
OY 


=") = Ts Pits 3¢ 


mevcwplied to the north south bouncary nodes for the terms 
Oeeet) in equeticn II1I-23. Integrating the inner preduct 
CB e°9 5.0, Dm Gpeirts pro@uces the Doundary terms 


2 : / . , 
jf v (BY 5)¥, dxdy = JP oo sts)¥, Coca 


I [v-(¥, v(8,7,)) = ADV .)> vv 4] dxdy 


7 a j ° 2 x A % 
§ ¥,7(9 575) dr I vV, v(9,7,) xdy 


V Piao 


ee AIG 4? * SVayel yy?! 


where zs is a unit vector normal to the domain and dr is the 


differential e¢istance along the path of integration on the 


perimeter cf the domain. 
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mee ene SoOntour integral in equation III-31 and put into 
salerkin form. in the sare way as in the one-dirensionel 
Beco vive equation in Chapter II. The resulting term is 


Gerivec as follows 
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f Virg.v Joa dr = $ ay v, ax 
Ax. 
= ( a = a5 
: OD sy + 2u, ra] be III-32 


Peneticn [11-32 appears twice in the continuity equation 
Tyr 22, fer time levels (n+i) and (1-1). AL1 values of wu are 
Known “or time (n-1), since they are saved from the previous 
Salcubations. Fowever, u(a+l) has not vdeen computed. To 
Smee fOr uln+l), doth W(n+i) and X(nt+1i) are needed. YV(ntt} 
meeselwea first from the vorticity eauatioa. X(n+1) needs 
Mamet) as mart of its solution and JZ(nt1) needs ulnt+tl) in 
foe “Solution. to avoid this oprobler, it is assumed that 
Minti) has a negligible contribution to the solution of 


u(z+i) and only Y(n+1) is used. 


oi 





IV. CCOMRUMOR TPPEEYeONT ATION 


RS cormulation anc general theory of the finite element 
method was presented in the previous chapters. The objective 
Meee S «€6Chepter is to discuss sore important comoutational 
SemectS LerFteining to the imrplerentation of the finite 
Semen. prediction syster. 

The main advantage that the finite element method has 
Omer other prediction techniques is its generality. 
Conceptvally, it seems possible by using many elements, to 
approximate virtually any surface with complex bcundaries 


* 


Me initial conditions to such a d€gree that an accurate 


ray 


meemebeon cada be oftaeined. In practice, however, obvious 
Smee neering limitations arise, amost important one being 
meeecest of the computation. As the number of elements is 
Biereeased, a larzer arount of computer time is requirec for 
mr omecast. furthermore, the limitations of the program and 
Pee computer may prevent the use of a lerge number of 
elerents. These limitations may be due to the corputer speed 
mime, storage aveilability, or round-off errors propagated in 
m= COmpitations because of finite precision arithmetic. 
miso, tem malfunction of a hardware component, if tne 
Damemctse) iS Carriec out using fWany computer hours to 
emennem™, can oe a serious proolem. It is therefore desiradle 


FO use efficient finite element programs. 
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he effectiveness of a program depends essentially or 


Mme es OllOWineg fectcrs. Firstly, the use cf efficient finite 


‘D 


iD 
da 


Ment tecnnrPeves is Maperta nt. 
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M@meero=mmidgg methocs eénd sophisticated use of the available 
memouter hardware and software are importent. The third very 
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Sclyememe Vceveloprent cf a finite elerent 


'cS 


proerem is the use of appropriete aumerical technicues. 

mee worticity-divergence model cescrited in the previous 
meepeer’ ic inrplemented on tre 13M 3223 computer located at 
the Naval Postgraduate Scheol. Some notatle features of its 
eeeckitecture are the three trillion bytes of virtual mass 
Storage, of which eight mega bytes are available to each 
mor, 2ad the 57 nanosecond macnine cycle tire. The model is 
executed rostly using a l2xi2 element domain requiring 42@k 
memes set Storage amd £9 seconds of CFU tire to execute. 
Bxceeding execution time and/or available storage is not a 
Seer l@m, im fact the system allowed a lot of flexibdility 
eupine the implementation pmase.of the rodel. 

Poe seurse code is written wsing FORTRAN IV and compiled 
meee Outimizing Fortran & compiler. Apvendix A contains the 
Semree cede listing. which is divided into five subdivisions 


@Peleneatin 


a) 


tae Boeiveal structure of the pregram. 


ma ricanaMy &RCEIPECTURs 


Ur 
y 
Oo 
JQ 
“a 
fu 
<3 
tty 
(D 
fu 
ct 
6% 
4 


es incorporated in the model are: 
1) Modularity. With only a few exceptions, ¢ach 


Boe ic jirited to one page of FORTRAN code. This makes it 





eesier to comprenend the program. Zach module performs only 


eee tesk -rCr exerple, subroutine CONTS0 comoutes the value 
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fee Pee oeine Verrer for tHe sontinuity ecauation. tLixewise 
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(sO One ™redule for the diverzence and vorticity 
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oes. SO itplement < new set of equations, only these 
mea:les would heve to de altered. 

Pee cs it ly comtroilabite switches. Switches ray be set 
moerueer print. plot or tabulate harmonic analysis data for 
[ee 6f he available fields. The ability to disvlay 
Paewermediete results allows each portion of the algorithm to 
te ronitored for computational adciustments. This also makes 
it easier for unfamiliar users to decome acauinated with the 
computational model. 

3) Forcing term subroutines. ic Devi ons 
implementations. each forcing term was calculated ty a 
Specie] subroutine. In this implementation, the calculation 
ere accomplished by general purpose routines which sirplify 
eae Peplementation.of the complex orosnost.iccequationms._This-. “<= 
allows implementation of different equation sets (i.e. 
Z2roclinic “odel) over tne sare domain with rinimal effort. 

4) Decufentation. Zach variatle is defined ty a 
ellert pWr2se (Appendix 4, A.). The function of each mredule 
ms eeliscriped in an introductory paragraph. Shared data is 


meecee tn tamece common blocks and identified with each 


cutroutire which uses them. A Sudroutine index is given. 








Be. Bein Prograr 

peel emain program is short, calling only six modules 
mumen ememwect the basic sequential flow of the wrodel. It 
Seerts with initialization of all model parameters (i.e. 
Meeel opntions, domain, finite element arrays, inner 
meecicts). ‘t then initializes the input fields (i.e. 
eecnctential heights, stream fuact Lon ase velocity 
Deamerneti al ) and iSueeoliemede by initialization of all 
réreining cepéndent variables. At this time the model is 
motally initialized and time integration begins. As 
mentioned previously two techniques are employed for time 
Bovesreation, eacn having its own module. Upon completion of 
mee last forecast, the program terminates. 

arrays are the only data structures used and are 
eroupec using iS different common clocks. Several arrays are 
Mee aS 6 Static €[CUCUClimk)6lListsSy 6as described in detail later, 


meee Simplified the algorithms. The common block forrat has 


t 


mm -crertasge ec recaucing thecoveral beenecution-time ofa the 
mueram., Fost of the arguments passed during a call to a 
Suoroutine are conrtainec in conmon. This requires less time 
Bemmemecute sijace no parameter passing is required for the 
Smeets, Another benefit of this rormat is that the coce 


vecores less cumoersome and more readable. Each variable and 


Pee 15 2efined in the first subsection of Appendix A along 
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Wien 2 Daze index for the subroutines. 


en 
Cn 





eee ielizationa Phase 
SMememcix Ae Section € contains all the Subroutines 
meee onring the initializeticn pnase of the ovrogram. Fror 
Mme user ooint of view, the most important suoroutine is 


MT -G 


tad 


, tne first subroutine called, which contains all the 


eye 1 variables Pant Somurol tee Girrerert -optioas 


(Wy 


maieabie per run. this is the only subroutine that is 


fw 


G@meesec tc run the different experiments, assuring that no 
Bea Computational technique is introduced. The selection of 
Getions are: 

1) - channel location - the channel may be 
smht ved nortn or south by onresetting the 
HOmmussoupn latitude limits in INITG3. 

2) - variable geometry - the node positions may 
te @rouned for more dense node patterns to 
yield higher resolution. Two variatles R1 
anc Re set toe ratio used to vary the 
spacing = alone tae Xeeweand= wey .° .axe $4 
respectively. 

3) -— imitial field wave length and amplitude can 
Beraiterec to produce varions effects. 

4) enmatge the initial mean flow. 

&) - é¢iffusion can te entered for any of the 
Taree progzostic equdtio1s. 

&) -~ maximum length of forecast period may 0»be 


saommen dae a Orint, plot or &arrcnic 
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analysis of any cependent variable may be 


requested for aay time interval. 


emcee wine experiment is determined, the options 
Listed ebove are set. The prograr is reacy to be executed. 

snrew largest part of the avy ati Za t Lon pnase 
memeests Of establisning the domain and producing all the 
Bares Clerent computational vectors that remain constant 
throughout the experiment. 

The first several steps in setting up the dorain are 
@emcoern@c With indexing. Subroutine CORRES is called “irst, 
where @ll the nodes (grid points) and elements (triangular 
areas) are numbered consecutively starting at the southwest 
Cormer of the domain and moving eastward across each row or 
mero vee. Yor “eath elementy- a record of all of its nodes 
Beeetices) are stored in array SLMENT (M,Z), where Mis the 
oeere! number of elements. To facilitate the inner vroduct 
Pveluetion later, a local numbering system is required for 
ech element. Thet is, fer each element, its nodes are 
feor@e CoUntérclockwise in a positive sense. The first node 
mewewe>, is arbitrary. 

Sitnm the domain divided and numbered, a cornectivity 
list (the corresvordence detween each node and the neighbor 
moees) ts constructed for each node by subroutine CORREL. 
B@en ode is adjacent to four or six other nodes depending 
@e wether it {s a boundery or interior atode, resvectively. 


These adiacent nodes, olus itself, make up the connectivity 
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Meee for O2e node. The connectivity lists are thea 
fommeaeeteved s€quentially starting with the first nodes 
Momeectivity iist into the vector NAME (NN), where NN is the 
mem time modes in each ccnnectivity list. (i.e. for a 
té2xle dorain with 186 nodes, and equilateral elements NN = 
ia42;. ‘or the first tire during the initialization phase, 
special ete nvion Po eewve4. “to mcyclic® continuity. 4s 
meeremose@ Parbiver, cyclic continuity is the joining of the 
Mest anc west boundaries to create a cylindrical channel. 
Spe comnmectivity list for these east/west bcundary nedes 
meer «|6lDeE)6C COMpleEte «6©6tol)6€O6hdimnsure§6proper continuity for the 
Sermombations later. 

"he connectivity vector NAME is frequently used 
Guri2 2 most computations. Twe utility vectors ISTART 
Weerrdaaime the starting locatica in NAMB for a particular 
node) and NUM (containing the number of nodes in its 
ceommectivity list) are used to locate and index through the 
mector NAME, ac will besseen shortly. This same: techaique.is 
moec to index taorough the coefficient matrices and used 
ime Mest of the node interaction computations. 

™ne ophysical properties of the CAaine! are 
mepcibeted ext im subroutine CSANAL. Eere the north and 
South latitude boundaries, which were pre-set in INITGE by 
mee “ser, are used to compute the grid spacing along the x 
Mme @8ls. Since this channel simulates a belt around the 


earth. the ragnitudes of both DEFLTAX and [ELTAY (mreters) are 





meeporteons! to the width of the channel divized oy tne 
memeeer of grid points in the y axis. 

The Cartesian coordinates for each node are comouted 
Pe sudrcutine LCCATE using the DELTAX and DELTAY calculated 
Py CeaNAL. If the option to use varying grid geometry is 
cesired, subroutine TRANS transforms the erid geometry. 
PRANS als? computes the corresponding new Cartesian 
meomarvete vwelues for each grid point and calculates the 
meotmur DELTAX and fDSLTAY within the domain. When the 
meemetry is changed to create a smaller DELTAX cr DELTA’, 
ome tWo cdimensitonal advective stability criteria is also 
meemecc. Anew time step DT has to be comouted using 
Seuvation ‘i1-2@. Since TPANS transforms the geometry, it 
eae> CoOmMputes the new DT. 

Seo wier treactormation is required as discussed in 
Eee ver TI. “he transformation from Cartesian coordinates to 
mua “Ccoorainhates ie meeded to perform the area integration 
o> the inner oproducts. Subroutine AREA computes these 
Preasformetions ‘exactly es outlined in Chanter II, Section 
meeeeenn Syclic continuity is very important and svecial 
oeame 16 meedec to insure’ oroper transfcrmaticn. 

eollowide the area transformation is the computation 
of é1] the inner products thet are required to solve the 
Cquetions. The advantage of using area coordinates is that 
the inner prodrets (function of space coordinates only) are 


Somemeea and stored cnce and used repeatedly without 





recalenlation. Subroutine INNER computes and stores these 
Meeeuets Using the formulas derived in Chapter II. 

“he coefficient matrix, dimensioned NxN, where N is 
tme totel number of nodes, is a matrix of coefficients whose 
pews ere the equetions of the system to be solved. As 
meeeussec in the computational téchnique section, the 
memcers Of this sperse fatrix are the inner products for the 
ert mend sides of the equations. Three coefficient matrices 
ome use@d in the solution of tne equations. The diagnostic 
Memetiogs ([LII-7 throwgh III-1&) wse a coefficient matrixz 
mmo 6 8DDelChC6Ucilcnmer©6product§3 «|6<V pons Pevbich is constructed by 
SerpPeutine AMTRX1 and stored in compacted form in vector 
Giwn) by subroutine ASEM3L. Eowever, when solving the 
mmoenoOstic equations. trese coefficient matrices have a Df 


‘time step in seconds) term, so that these matrices are not 
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Smebed until the time integration begins. The vorticity 


fu 


fu 


me Givergenece equations (I1I-24,111-25) use the coefficient 
Beerin =(0N) with tener uproducts V gy Vay 2, + VayWy? in. 
Sou viae the -oiscson equations for the stream function and 


Bebocity potential, respectively. The continuity equation 
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emesis «6C8lCUCCOMbEaation of ianer products in its 


coefficient matrix F(NN) as fZollows 
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Be solve the Felmholtz equation. At the start of each time 
integsretion rodule, Svtroutine AMTRX2 is called to construct 
mee temo Coefficient ratrices © and F. 

shese tanded and sparse matrices are compacted into 
Meo“ OrS tO save storege during their assemblage by ASEM3L. 
The vectors are dirensioned NN, as is NAME, the connectivity 
meeewor, a2c both use ISTART and NUM to index through ther. 
This compaction routine was used by Kelley (1976) and Older 
‘19Si) in théir models, but was developed by Hinsman (1975). 

To illustrate matrix asserblage using an element Ddy 
elerent technique, consider Figure 8. Note that this 
illustration is for element number 5, but all elements are 
Bmeeted 12 @ similar fManer. The compnutational technique 
required that for each peint (node) describing element 3, 
mamMely nodes 2, 3 and 14 stored in array ELMENT, the inner 
Pod yc t ee? between those points te distributed to their 
opemer boication in the coefficient matrix. 

Subroutine AMTRXD-batlds the inner product wodalz 
Beeemection and stores it in matrix 3, dimensioned 5x3. 
Fie:re & illustrates the 3 matrix for elerent 3, where the 
inser prexdwe t <V 5 no > is the multiplicand of the 
eewresvpomadle@s besis ard test functions. respectively. 

Beeeeeocah dispensing of interactions is doae in 
ASE“EL. Consider the second row of (3] in Figure 8. These 


are the interactions between node 3 of element 3 to the test 





wemex - builds (3] ror ore element, passes (3] and the element number 
associated with (2] to ASEMSL. This process continues till all element 


node interactions are assembled in the coefficient natrix. 


Yovo VoVy (VoVan 
B = Vay Way, Vo ELMNBR = 
se 5% 5% 3 
Vouo Vants  “se¥ry 
AScMBL - assembles the node interactions into the coefficiert matrix 


(c] » which has the same structure as MAME. The following diagram assembles 
inner product Yaa into the coefficient matrix (c], for element 3. 


AME 






a 
a 
287 


288 sue |133 1145 





155} 1035 | 
156 |1040_ 10 
ASEMBL pseudo-code START(3) +11 
do 1 = 1 - 3 z 
44 = BLMENT(L) ( 3 t3 
14 
do j = START(11) -—> LAST(i1) LAST(3) #15 


jj = NAME(j) 
do k =1l1 3 
kk = SLMENT(k) ( 3 
14 
at He = Gg ) ~ then 20 
C(j) = c(3) » B{4,%) 
else continue 
pe do 
end do 1043 
end do 104d, 





Figure 8. Assemoling and storing the coefficient matrix for element 3. 
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functions. ASEMBL locates nodes 2°5 connectivity list in 
Seeeeeciee ISTART and NUM. In Figure 8, this list is 
delineated by START(3) and LAST(2). Now ASEMBL steps through 
meereoounectivity list for three iterations. During each 
pess. ASs”2L ts searching for one of the three node numbers 
fer Clement ¢. When a match is found with one of element 3’%5 
memes ‘i.e. 2,3 er 14) and mod@e 3°S connectivity list (i.e. 
fe, 14,15 or Ali Vene orieper position, to which this 
interaction is to be added has deen found in the coefficient 
Metrix. Since AME and vector C, the compacted coefficient 
matrix. are dimensioned identically, the same pointer (i.e. 
fee2 Piewre €@) is used to index through both arrays. This 
Peocedure is repeated for ail elements in the domain to 
eememole the coefficient matrix of the equations. The 
Pomewen-code for ASBMBL is sneown in Figure & to facilitate 
SeeGppine thrower this example. 


The domain and all finite element work vectors are 


me the ligicgeae this moint. Soapbroutime ERMSST. is called: later. 2: 


meeecommaa.e interpolation voints for the harmonic analysis 
cmeromtines. 

Mee Pest bhase of the initialization process 1s the 
imitialization of the dependent variables. The tnree input 
Pields geopotential heights, stream function and velocity 
potential are computed in subroutine IC using the equation 
Peet ii-oG. SoyPver, the variables calculated from the 


@ia2znostic equations neve to be computed using the input 
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These variables are used during each tire step while 
Eeamome the vrognostic equations. 


nm 


Meee Cisgnostic eqwations are solved in subroutine 
feewer, first during the initialization phase and later 
@erees the time integration phase. Fech diagnostic equation 
calis its own module to compute the value of the forcing 
eee tion anc stores the computed valves in the vector RES. 
Thesé eauations all use the same coefficient matrix when 
Moeving the diéenostic equations. Subroutine SOLVER is 
sufficiently genereal to solve each equation. SOLVER uses 
Wector 225 and coefficient matrix G to under-relax towards 
mee solution. As mentioned previously, the coefficient 
matrix is strongly diagonally dominant so that three passes 
ewer the domain are sufficient. At the end of DEFVAR, output 
me generated depending on what print, plet, or harronic 
analysis controls were preset. 

This corpletes the initialization pnase of the model 


enc the program for the forecast phase will be described 


&. Ferecast Phase 
The “forecast chase is accomplished in two steps. The 
first time step is made using two half steps by subroutine 
mereno. H@Fe the vorogmostic coerficienrt Matrices are 
constructed using half the DT value dy calling AMTRX2. 
BMTEIe wees the same computational technique to construct 


tne coefficient matrices as described for AMTRZ1. 
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beweh of the prognostic eauations III-23, III-24 ane 


Meeoegecaiis its own subroutine (CONTEC, VORTFO and LIVEC 


2 


espectively' to compute all the terms on the right hand 


in 


mee, eich are stored in the vector RHS. After computations 


y 


ty 
oO 


me 


> are completed, subroutine RELAX solves the equations 


> 


me cwer-relexaticon as described in the computational 
Deemmpane in Chavter 3. Once the solutions for the (n+1) 
meme, SUED are completed, DEPAR is called to update the 
Variables from the diagnostic equations. Twe passes through 
MATZNO advances the solution fields one time sten. 

The remainder of the forecast period is integrated 
Meeme the leapfrog scheme. Sudroutine LZAPFR performs this 
matveeratvios using the identical format as MATZNO, except 
Meee 0 equals two DT. At oreset times as specified in 
Mer@ee, the different fields are saved for printing. This 


process continues until the final forecast time is reached. 


tat 


UFILITY PoOrlrubes 

Bee (me @audtion formulation is completed, as in 
eeomseer ‘IT. all the ianer products and types Ot 
Seeeerdtions are known. Versatile modules can be written to 
meeoomr these computations. Consider a term of general form 
“a v; Vs? where i is the node about which the term is 
evalueted and the j’s are the nodes cennected to node i, or 
fee SUrrounciage aocdes. The inner product values Sen are 


@iready computed end stored for ell the nodes, during the 


Pemeiaiization chase of the rodel. 


O) 
cn 





The reméining computation to complete tne evaluaticn of 


Dais t€rr is the multinlication of the scalar coefficient of 


£ 


fu 


femeoce i wita the corresponding innér product Sin Nae Stone: 
node ‘. “his requires indexing through node i’s connectivity 
fer oucre@ in vector NAME, and for each node in the list 
Mme iviy ane total the products. The cumulative sum of these 
Teltiplications is assigned as node i’s contribdution for 
this term. Subroutine TER”3 performs this exact computation. 
mit thet is passed to TERMS is the scalar field A and the 
sien o? Mee faner coreduct, TERMS then computes the 
eemeritution for each node in the domain and accumulates it 
mo tHe work véctor 25S. 

snree other utility modules are; TERM1, which computes 
Pee First Scalar multiplication for triple inner products 
(die. “A; ¥53,V,,V,>). The product CVV WD, is again already 
corduted and stored by subroutine INNER. TERMI computes ae 
Mae 5 > ard constructs a compacted vector sirilar to the 
coe ficient matrices. This weduces the effort of multiplying 
the second scalar to a TERMS computation. TERM2 computes 
meee interaction of the following type <A Vix Vix? wnere both 
the basis and the test functions are cerivatives. Lastly, 
Suvroutine TERM4 computes node interaction for terms as ce 


See aereme only the basis function is a derivetive. 


When examining the right hand side of the equation sets 
memeee, ifi-2f and II'-25, it is obvious this implementation 


me a subscripting nightmare; however, the use of tne utility 


E6 





modules TERML. TERM2, TERMS and TERM4 simplifies the 
igplementetion to only determining what order to call the 
meee ity wedules. ExXamination of subdroutines CONTEC, VCRTEQ 
mae TIVWE, which compute the right hand sides for III-23, 
Pil-2@ and IITIi-25 respectively, illustretes this fact. No 
other subroutines or caiken) at pons were reauirec. 


Itplementation of these equations required minimal effort. 





No See Live MOlaL EXPERIMENT 


The orevious tWoe Chapters presented tke detailed 
Sormubctior anc implementation of the vorticity-divergence, 
Eeeelow water equation model. The results from this rodel 
rl eoe Comperec with tne results from the primitive model 
in Gheveter iim Ordotritete miteroretetion of the 
comparisons, a brief descrintion of the primitive mcdel 
Pouleaws. oee Selley (197€) for a detailed discussion of the 
Satire wodel. 

Also DRe@sented in this chapter is an experiment which 
@emonstrates sienificaat improvement of the solution from 
the primitive model. Kelley’s implementation used elements 
which were right triangles. Clder (i¢$81) showed that 
equilateral elements are far superior Pe triangular 
elerents. - This experiment re implements the vrimritive model 
MSsing@ @quilateral elements and a comparison 1S Mace between 


eee results of both implementations. 


mm. *ClEL DESCRIPTION 

mueeorme of the csarctropic, shallcw water, primitive 
equations developed by Phillips (1959) is used as the 
Rovernineg ecuétion Set tor this medel. In UCUartesiar 


Seoermcingates tie Cavation set 1s 


eS 





Og 3 3 

se se Oe) - ee (Ve) V-1 
a Ox ele 

du ag au du 

aaa “Ae oe y — sey, Vee 
at x Ox oy 

av ag av QV 

—. ee ees Bl ee Vv a et. Y-3 
Oe dy dx Oy 


she finite element formulation of this set of equatioas 
evaluates the height and velocity components at the same 
Meee peimts. This is an important consideration, because 
the other models (the linearized model, see Chapter VI, and 
the vorticity-divergence model, see Chapter III) either 
Stegeer the cependent variables or have the property of a 
Stazgzered formulation. When comparisons are made between the 
moeemrs. it is this lattice structure that is being compared. 

wees eormm™ 6 6CtlhCUtme 6€6CScmalicw ©6 water equations includes 
gravity waves as a solution. Gravity waves have a maximum 
phase speed of about 220 meters/second. When the correct 
fime sten is eelculated using eoimeicwon  Lll-26,5 4 
pommeacerebly <faller time step is obtained compered to the 
Peo tame step permitted in the vorticity-divergence 
Memembhetion. this is an important feature. If solutions from 
ell models are equally as good, the dest formulation would 
me cec@rrined using the computational time required to 


Procuce the decired forecast. 





All mocels use the same domain structure. In fact, tha 
meer «6 eCSCTibed)«€6in)6€6(Chapter)6«€6[1I)6was patterned after the 
memein implemented dy Kelley. Agzain, this domain simulates a 
Meee eroume the earth, with cyclic Coat las by Which 
memanawes the east-west boundaries. Rigid boundary walls 
Beest ct the equator and at 39 degrees north latitude. The 
memamn 1S composed of a 12x12 point mesh and subdivided into 
the right triangular elements illustrated in Figure $. 
Rotice that the grid points are not shifted as in the 
Eau@ilateral element implerentation shown in Figure 6é. 

7he following boundary conditions are imvosed: 

1) = 20 Snoisis  ‘Chamael flow at the latitude 
Soundaries. 

2) - @ geéeostrophic balance at the channel walls 
Imoosed on Wire continuity equation V-l. 

mis Podel has @ sifplie second order diffusive term in 
Mee Gauatiois of motion V-< and V-S. However, for the 
Durpese~-of e€valueting these different-forrulatfors, tnits 
mereroa Wes 20ot implemented duriag the comperison phase. 

Meets conditions consist of a singlé wave in the x 
memeetiond did a half wave in the y direction. The initial 
Fi€lds "or the three dependent variables are shown in figure 
foe fee Wexisum zonal wind perturbation of 5.5 meters/second 


is superimposed on @ mean zonal flow of 12 meters/second. 
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wees 10. Initial fields for the primitive model. Both the x and y 
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axes are multiplied by 10° Km. 
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21S cursory description of the primitive rodel mentions 
fee DOSS Sienificant features that will weigh heavily in 
mee Comperisions later. The Galerkin implementation of this 
mocel 1S Similar to that presented in Chapters II and III 
anc) 6€6©6the)6system of ecuations is solved using a Gauss-Seidel 
mmerative procedure. ‘Further details eo1cerarmie this 


primitive model are given in Kelley (1976). 


e. ‘RESULTS 

“Mis @xperimert involves the shape of the slements. As 
mentioned previously, Xelley’s implementation subddivides the 
momeen into right triangular elements. as tllustratec in 
Beeure S$. Considerable smali-scale noise was observed by 
merley im the 48 hour forecast solution. 

me Gransition from right triangles tc equilateral 
eremreleés chanses the size of tne domain. With right 
triangles. the«sx anday grid spacings are equal (380 KM). A 
meree grid matrix has a length and width of 3600 <*> With 
Ser, lateral triangles, the+zr and sy grid spacings are 10 
meereer “equal. Arbitrarily, the ay grid spacing is hele 


constant (390 £M) and a newax grid spacing computed by 


ax = ay/cos(32) V~4 

A 12x12 grid matriz with equilateral elerents nas a width of 
eee $f anc a lezeth of 4245 KM. 

Figure 11 contains the 48 hour forecasts produced using 


both types of elements. The three dependent variables fields 
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nd equilateral triangles for element shapes. 


amplitude of 


~= 
a 


48 hour forecast comparison from the primitive model using 


2 
both right triangles 


igure 


w 


> 


arbitrary perturbation 





mee Compareé for each. The small-scale noise that felley 
ebservec is present. The three fielés shew varying degrees 


~*~ 
> 


mearcistortions, which are especially noticeable alorg the 
Poundéeries. 

Clder (1981) found that theroot Tean Square error (RMSE) 
wes reduced 22 percent by using eauilateral shaped elerents. 
“his ifvrovement is apparent on viewing Figures 115, d and 
#. ‘the contours are smooth and the dvcundaries are 
noise free. Felley showed excellent treatment of wave 
eroperetion by this primitive model. The lowest resolution 
grid ‘'€x6) tested by Xelley was within four percent of the 
true phase velocity. Changing the element shape nad no 
apparent effect on the phase velocities. 

Fecause the outcore cf this experiment was a 
meee t icantly improved forecast solution, future comparisons 
Mece the primitive model will be made using equilateral 


elerents. 
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a) Piecewise linear basis function. 
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b) Piecewise constant basis function. 


Figure 12. 


Different shaped basis functions. 





BR. <:CUATION PEFORMULATION 

-he primitive forr of the shallew water equations 
oreséented in Chapter V (V-1,V-2 and Y¥-3) is used to derive 
the lineérized equations needed for this experiment. The 
memoecity Catiations V-2 and V-3 rerain unchanged and a linear 
Mesis function (7 ) is used to epproximate the u and v 


verieables. 


ere comtinuity equation V-1 is linearized as follows: 


a2 af 3j gus 
—=- -vp -9( —+—) VI-1 
ot z y Ox oy 


negligible 


where 4 is the average gzeonoteatial over the domain. A 
viecewise constant basis tTunction (Ws iets used mmio 
@cproximete the geopotential. aes linearization is 
reasonable in this case because the Rossdy radius of 
@eforretion be is much larger than ax [see Williams and 
Zieckiewicz (1981)].-The. Galerkin-method-is applied to. this 
lineérized equation set using a piecewise linear test 
memction fer 4-2 and V-3, and a piecewise constant test 
Sumiction for VI-1. 

Boesoiecemise comstamt basis function has the proverty 
Oe» €isplacing the geopotential to the centroid location of 
the elements, which should give the same effect as 
Steegering the grid points, as seen in Figure 13. The 


Sensity of geopotential data in the dorain is now greater 











o 
* 
= 
_ 


Mo ¥ u,v U,V u,v 
Cc °C — J 
é = . 
g M) g g 
u,v u,v u,v u,v 
6 3 


~~ 
iy 
fa 
on 
ee 
ar 





u,v u,v u,v u,v 
oD = -- CO Oo 
- D g g 

‘ 
u,v Uv U,V u,v 
gO 





Figure 13. Staggering of the geopotential 
about the velocity. The o symbol are the 
actual grid point locations and the ® symbol 


the centroid location of each element. 
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Beer XRE censity of the velocity data. Instead cf a 
femererenk tal value for eack node, there is one averazed 


Each elerent. 
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eemeaaet form of the Galerkin equation for VI-t, vV-2 


eee v=o atter verforming the time differencing is: 


eye oe 7 eerie ni n Names: ae 
a <9 .;#4> = @ XH 56 W,> 4td iu er “wy ee \ ia 


STi 
n+l n-l n noon 
Uu @ <@ a s ; s We ar f Cw V = e a V 
3 b V;> u SE at ld, cs i? DTV Vi VG 
od ee -— j 5 —_ ae Wipes 
ee cy We> = wia ee - We ee. SS + Cy. ov... 
Dili el ee se J kk A 
mn nN 


Pais !ipearized set of equations VI-2, VI-3 and VI-4 is 
Semevec Using a sauss Seidel iterative orocedure. It is worth 


Menvionir2s thet the coefficient matrix <W a oy? in equation 


~} 


- 
t 


Bemis agile ston zero coefficients equal to one, since the 
Breese retion of the inner product C5 Wy > involves piecewise 
@eactant functions. 


ans eouation set is implemented rather than the 


b ty 


Sowvetion set '-1, ¥Y-2 and V¥-3 using ell oO tne existing 
Dmetawe 3 moce}] code. The major modification involved the 
Way Unet the geopotential was referenced. wWith an average 
Meena teat over each element insteac cf a value at Sach 


meee for ¢ 12x12 Wesh, there «re 288 Heonotentiel oacinrnts 


Weesws tee 1L5& velocity (u,v) points. 


Q) 
ee) 





t+ 


=. S2cULTS 

Pere results from this linearized model are compared to 
ee fromr the primitive model. The initial field fcr this 
experiment, rigure 14 has a maximium perturbation zonal wind 
met is one fifth of tne value used in Chapter V, Figure 14. 
The Tear zonal wind remains 14 meters/second. The 48 hour 
Porecast soluticns for each field are compared in Figure 15. 

feat alternative formulation shows some promise, 
although tere are some minor perturbations in these 
Somtours comparec to those in the primitive solution. No 
mmelenation is offered for this small-scale i1o0ise, although 


meecibly the other formulations presented by Williams aad 


Zienkiewicz (1981) weuld improve the soluticns. 
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Sigure 14. Initial fields for both the primitive model and the 


linearized model. 
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Figure 15. 48 hour comparison between the primitive model and the 


Linearized model. (APA = 1.1 m/s) 
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me. MOET ICITY-L1VaRGENCE MCDSL EXPSRIMENT 


ms) 23 the vreéevicus two chapters, a comparison tetween 
two mocels will be ovresented. The results rrom tae 


Meeticity-civerzgence rodel are comoared to the results from 


9 
e7n 
wd ah 


(D 


Peemoweewe motel, i909 fully exploit the differences 
between the mocels verformances and indicate the streneths 
ae2 weaknesses or eacn formulation, three domain geometries 
mae 6 6tUnr ee $Aitieal conditions are sed. All Solutions are at 
Memeo >. xcept for one case which was extended to 96 hours. 

eon thiese two fiaite elerent Poumu latio7s core 
Meesiond! insight is osdtainec concerning tne execution time 
meeurred as tne erid resolution cnenges. Lastly, a brief 


MmeorrcstioOn on the sce@msitivity of the computational technique 


me eGo EOMAINS ANT INITIAL CONDITIONS 

The three domain geometries used in the model évaluaticna 
Meme ies rated 12 Fireure 1. All domains consist of a lexle 
elerent resh with e2quileteral shaped elements (15€ grid 
feepic COmememtrt, is imaosec on tae eest and 
west Doundaries. The domain hes dimensions of 4245 AM along 
mem xeenmic anc 2575 4M alone the y axis. 

Tee Rmeieuler domain (Figure Gay “Rds o) ea nS ce 


Seeerabutt!on of @rid voiats, with a mMminimur grid spacing 
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Figure 16. Test domain geometries. 
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figure 24, Initial fields for both the primitive and the vorticity- 


divergence models using the smooth domain. There are two waves embeded 


in the flow and APA = 1.1 m/s. 
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Figure 25. 48 hour forecast comparison between the primitive model 
and the vorticity-divergence model using the smooth domain. Two waves 


anee Smeeced in the flow and APA = 1.1 m/s. 
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Figure 26. Initial fields for both the primitive and vorticity- 


divergence models using the abrupt domain and perturbation amplitude 


of 141 m/s. Note that the 


element shapes are not equilateral triangles. 
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Figure 27. 48 hour forecast comparison between the primitive model 


and the vorticity-divergence model using the abrupt domain. 
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Comparison of the computaional times for a 48 hour forecast 
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a octry aT TONAL SENSITIVITY 

me= Section will offer an exvlanation for the skewing 
oe Contours in the SE huur forecast solution over the 
aeootwe deWein ‘8ieure 23;. Aas mentioneé ovreviously, there 
wes met ©Pnough Lime ewaileble fer “Fine tunire the 
m_eemmeramation coefficient, which is used while solving the 
system of equations. The overrelaxation coefficient may de 
Bere soensitive anid small changes can, on occasion, radically 
M@eenee the rate of convergence. The opntimal value or the 
Meerrel=x%eti0n ccefficient depends on the specific ferr of 
Mee e@erriciemts of the equation and the error distritutiom. 

P= eEcuation set to ve solved consists of three 
Pametaonrs and each equation reauired its owa relaxation 
ee"ficient. When solving the ecuations over the resule 


Germain, the entire system is well tehaved and an optimur 
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Memeretion Coefficient is easily determined. Eowsver, as the 
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Seeger sine tuning. 
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Bmore nour forecast uses tne smooth domain. The 
Memeo retitraeineae] eszrid points are cempacttecd, creates a ie®rser 
ema TNE ritddle sf the channel. The eceefficients 
erizinelly computed using the See ae domain needed 
meejustmenats to vrererly solve the equations. 

Demecimeyrate the significance for fine tuning the 
Mememmrron coefficient, consider the series cf plots ina 
Bereure 22. “he vorticity divergence eauation set can be 
Semoliziec by assuming the flow is non-divergent, so that 
mew tre worticity eouation naeeds to be solved. figure ca 
memo se dour worticity field using this ecuation over the 
member “omein with an overrelaxation coefficient of 1.3. 
Beet ield i¢ well defined with Srooth contours. 
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Figure 28. Computational sensitivity using the 48 hour vorticity fields. 
Fig. 28 A) the 48 hour forecast using the regular domain, 28 B) the same 
forecast using the smooth domain relaxing in one direction, 28 C) same as 
28 B) except alternate the direction of relaxation after each pass over 
the domain, and 28 D) same as 28 C) except the relaxation coefficient was 


Sia] tuned foré 1.3 to 1.297. 
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mere tee relaxation coefficient from 1.2 to 1.267 


Meo Wees tse much imprevec soiution in Figure 2&d. If the 


in 


mmcomeomroGg  coetficien: POimetibemouner two ecuations could 
Pee- eee eyerre tunmec, it anpears that improvece solutions would 
result. There are also other relaxation techniques availeble 
meer neve poteatia: for improvine tine solution while also 


Somreoere 2S at a faster rate. Some of tnaese technicues will 


me vested in tre futuer using this rocel. 





This reséarch investigated different finite element 
PorM@@letions for rae Snallow water equations. The 
two-—cirensione!? domain was a channel which sirulated a telt 
Browee tie Garton. Anal#tic initial conditions were used to 
Semelwiy the comperisons. Two formuletions were examined; 
Omer wetste cgifferent shaped basis functions and the cther 
meere 2 €itfPerent form of the Equation. Fach was compared to 
ime “oremivive form of the shallow water equations that was 
develoved by Kelley (1975). 

The use cf equileteral shaved elements which was 
Sesestec by tr. MeJ.F. Cullen sianificantly improved the 
Soluvtons comeared to Xelley’s model, which originally used 
meet teiemel@s as vsasis functions. Most of the other 
stuciés in this thesis used the e€aquilateral trianeagleées. 

williams and Zienkxiewicz (1581) svgvested the use o*% 
piecewrse Linear basis functicns for the velocity field and 
Beeeewr See comctant muwnmetions for the heizht ?ielec. This 
formulation was tested with a linearized COm vl atic 
eduation. ‘the results were peorer tnan those obtained with 
Téelley“s model. 

Most of the effort a this thesis was devoted t9 


iMeeewrenttetamhe end testing 42 verticsity-diverzgence mecdcel 





peur tC the on@s develonec by Staniforth and Mitchell 
(1977) and Cullen and Hall (1979). Several tests were 
peeee@atec which compere this formulation with Selley’s 
jemel. it Was found thet this mcdel executes anproximately 
Ome Orcer Of Magnitude fatter than does the orimritive 
Powe lation wsed “sy eKeley. secondly. as the spatial 
Mesorution wetween sridc points decreases, this formulatioa 
Omeewees a solution that is far superior to the primitive 
Momr. © disadvantage is its computational sensitivity, which 
Meemoires fide tuniag in solvine the elliptic equations for 
certain peormetrmes. It also requires 25 percent fTore 
feomet.er Storage, due to the more complex equation set and 
Phe. 2@ditional ‘Variables that ere treated. 

Peopmementation of fidite element models is not easy. 
Pomever, there is a lot of generality and redundaacy 
imeeccea in the computations. VYersatile modules were written 
Peer shen cantivy reduced the effort in implementing tas 
vorticity-civerzence medel. 

Further research is suggested using this finite element 
Poemuletion. It has accurate onase vropagation, is abdlé to 
hancle variable arid geometry, recuce the small-scale roise 
See) @ecrecse the model s executicn tire. Svecifically more 
aaerebicee WeEtTI0ds of solving the elliptic eauations shouid 5e 
imvestigatecd. Finally, the formulation should te tested witn 
Smmell~<eale forcinae, where its advantages shoulcd be most 


evident. 
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